Computer Implemented Method for Modelizing a Nuclear Reactor Core and a Corresponding Computer Program Product

ABSTRACT

A computer implemented method for modelizing a nuclear reactor core, including the steps of: partitioning the core in cubes to constitute nodes of a grid for computer implemented calculation, calculating neutron flux by using an iterative solving procedure of at least one eigensystem, the components of an iterant of the eigensystem corresponding either to a neutron flux, to a neutron outcurrent or to a neutron incurrent, for a respective cube to be calculated. 
     A control parameter is varied to impact a neutron eigenvalue μ through a perturbed interface current equation and drive the neutron eigenvalue μ towards 1.

This application claims priority to European application EP 09305767.7,filed on Aug. 18, 2009, the entire disclosure of which is incorporatedby reference herein.

The invention relates to a method for modelizing the core of a nuclearreactor, especially for calculating neutron flux within the core.

BACKGROUND

The results of such a modelizing method can be used to prepare safetyanalysis reports before building and starting a reactor.

These results can also be useful for existing nuclear reactors andespecially for managing the nuclear fuel loaded therein. In particular,these results can be used to assess how the core design should evolve intime and decide of the positions of the fuel assemblies in the core,especially the positions of the fresh assemblies to be introduced in thecore.

Such modelizing methods are implemented by computers. To this end, thecore is partitioned in cubes, each cube constituting a node of a gridfor implementing a digital computation.

Usually the steady-state diffusion equation to be solved during such adigital computation amounts to:

$\begin{matrix}{{\underset{\underset{removal}{}}{\Sigma_{Rg}^{m}\varphi_{g}^{m}} = {{{\lambda\chi}_{g}{\sum\limits_{g^{\prime} = 1}^{G}\underset{\underset{production}{}}{{\nu\Sigma}_{{fg}^{\prime}}^{m}\varphi_{g^{\prime}}^{m}}}} + {\sum\limits_{g^{\prime} \neq g}^{G}\underset{\underset{inscatter}{}}{\Sigma_{{gg}^{\prime}}^{m}\varphi_{g^{\prime}}^{m}}} + {\sum\limits_{{u = x},y,z}{\frac{1}{a_{u}^{m}}{\theta_{gu}^{m}\left( \underset{\underset{incurrent}{}}{j_{gul}^{+ m} + j_{gur}^{- m}} \right)}}}}},} & (1) \\{\mspace{20mu} {{with}\mspace{14mu} \left\{ \begin{matrix}{\Sigma_{Rg}^{m} = {\Sigma_{ag}^{m} + {\sum\limits_{g^{\prime} \neq g}^{G}\Sigma_{g^{\prime}g}^{m}} + {\sum\limits_{{u = x},y,z}\frac{2c_{1{gu}}^{m}}{a_{u}^{m}}}}} \\{{\theta_{gu}^{m} = {1 - c_{2{gu}}^{m} - c_{3{gu}}^{m}}},}\end{matrix} \right.}} & (2)\end{matrix}$

where λ is a first neutron eigenvalue, m is a cube index, also callednodal index, G is the number of neutron energy groups and g, g′ areneutron energy group indexes, u is a Cartesian axis index of the cube,Σ_(ag) ^(m) represents macroscopic absorption cross-section for the cubem and the energy group g, Σ_(fg) ^(m) represents macroscopic fissioncross-section for the cube m and the energy group g′, Σ_(gg′) ^(m)represents macroscopic slowing down cross-section for the cube m and theenergy groups g, Φ_(g) ^(m) represent neutron fluxes, such that theΣ_(Rg) ^(m,)..Φ_(g) ^(m) represent the reaction rates for thecorresponding reactions (absorption, fission), ν is the number ofneutrons produced per fission, χ_(g) is the fraction of neutronsemerging from fission with neutron energy g, a_(u) ^(m) is the width ofcube m along Cartesian axis u, and

with the relationship between the neutron outcurrents j_(gul) ^(−m) andj_(gur) ^(+m), neutron fluxes Φ_(g) ^(m) and neutron incurrents j_(gul)^(+m) and j_(gur) ^(−m) defined by:

$\begin{matrix}\left\{ \begin{matrix}{j_{gul}^{- m} = {{c_{1{gu}}^{m}\varphi_{g}^{m}} + {c_{2{gu}}^{m}j_{gul}^{+ m}} + {c_{3{gu}}^{m}j_{gur}^{- m}}}} \\{j_{gur}^{+ m} = {{c_{1{gu}}^{m}\varphi_{g}^{m}} + {c_{3{gu}}^{m}j_{gul}^{+ m}} + {c_{2{gu}}^{m}j_{gur}^{- m}}}}\end{matrix} \right. & (3)\end{matrix}$

The coefficients c_(igu) ^(m), with i=1, 2, 3, are characteristic of thecube m and depend on nodal dimensions and macroscopic cross-sectionsΣ^(m).

FIG. 1 is a schematic representation in two dimensions of a cube mshowing the neutron incurrents j_(gul) ^(+m) and j_(gur) ^(−m) for u=x,y and z; the neutron outcurrents j_(gul) ^(−m) and j_(gur) ^(+m) foru=x, y and z; and the neutron fluxes Σ_(g) ^(m). Indexes l, respectivelyr, refers to each left interface surface, respectively each rightinterface surface, of the cube m for the respective Cartesian axis x, y.Indexes +, respectively −, represents the orientation from left toright, respectively from right to left, for the respective Cartesianaxis x, y.

The steady-state diffusion equation (1) is also named NEM equation, forNodal Expansion Method equation.

In the state of the art methods, most of the computational efforts areconcentrated in the part dedicated to the iterative solving of a largeeigensystem corresponding to the steady-state diffusion equation (1).

In order to lower these computational efforts and therefore acceleratethe solving of the eigensystem, Coarse Mesh Rebalancing (CMR) procedureshave been used. In these procedures, neutron fluxes and currents for agiven iteration are multiplied with a corrective factor before pursuingsubsequent computationally expensive iterations. The multiplicativecorrection serves to suppress the presence of a non fundamentalwavelength part of eigenspectrum with the first neutron eigenvalue λclose to an exact value λ_(exact).

However, the acceleration effect realized in this way depends on thenumerical proximity of the highest coarse mesh level in a multi-levelhierarchy to the full-core diffusion level. Such CMR procedures maytherefore lead to very slow convergence or even convergence stagnation,thus increasing the computational effort.

SUMMARY OF THE INVENTION

An object of the present invention is to solve the above-mentionedproblems by providing a nuclear reactor modelizing method which offers abetter convergence accuracy, a better computational robustness and abetter computational efficiency so that relevant neutron fluxcalculations can be obtained within a short computational time periodand with a very good convergence accuracy.

The present invention provides a computer implemented method formodelizing a nuclear reactor core, comprising the steps of: partitioningthe core in cubes (10) to constitute nodes of a grid (12) for computerimplemented calculation, calculating neutron flux by using an iterativesolving procedure of at least one eigensystem, the components of aniterant of the eigensystem corresponding either to a neutron flux, to aneutron outcurrent or to a neutron incurrent, for a respective cube (10)to be calculated, wherein a control parameter is varied to impact aneutron eigenvalue μ through a perturbed interface current equation anddrive the neutron eigenvalue μ towards l.

The present invention also provides a computer program product residingon a computer readable medium and comprising computer program means formiming on a computer implemented method.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will be better understood upon reading of the followingdescription, which is given solely by way of example and with referenceto the appended drawings, in which:

FIG. 1 is a classic representation of the relationship betweenmacroscopic cross-sections, neutron incurrents, neutron fluxes andneutron outcurrents, in a modelized nuclear core reactor,

FIG. 2 is a schematic view illustrating the partitioning of a nuclearcore, and the association of driving factors with interface surfaces ofa cube according to an aspect of the present invention,

FIG. 3 is a schematic representation of the decomposition of isotropicand anisotropic parts of neutron outcurrents leaving the cube of FIG. 2,

FIG. 4 is a schematic representation of a multi-level V-cycle from theiterative solving procedure according to another aspect of the presentinvention,

FIG. 5 is a detailed view of box V of FIG. 4,

FIG. 6 is a schematic view illustrating the partitioning of a nuclearcore according to another aspect of the invention, and

FIG. 7 is a set of convergence curves of eigensystem iterative solvingprocedures for the state of the art method and different aspects of themethod according to the present invention.

DETAILED DESCRIPTION

In the following description, the case of a pressurized water reactor(PWR) will be considered, but it should be kept in mind that the presentinvention applies to other types of nuclear reactors.

In a first step of the computer implemented modelizing method accordingto the invention, the core of the reactor is partitioned in cubes 10(shown on FIG. 2) as in the state of art methods. Each cube 10corresponds to a node of a grid or network 12 on which numericalcomputation will be implemented through the computer.

In order to ease the representation, the grid 12 is shown on FIG. 2 asbeing two-dimensional, but it should be kept in mind that the grid isactually three-dimensional in order to represent the whole core.

The neighbours of the cube 10 with the cube index m (in the center ofFIG. 2) are the cubes 10 with the respective cube index m′₁(m), m′₂(m),m′₃(m) and m′₄(m). In the following of the description, the cubes 10will be directly designated by their respective cube index.

In a second step of the computer implemented modelizing method accordingto the invention, the neutron fluxes Φ_(g) ^(m) within the core will becalculated by the solving of an eigensystem corresponding to thesteady-state diffusion equation (1). To this end, an iterative solvingprocedure is used.

A removal operator {circumflex over (R)}, an inscatter operator Ŝ, aproduction operator {circumflex over (F)}, and an incurrent operator Ĝare defined, from Equation (1), by:

$\begin{matrix}\left\{ \begin{matrix}{\left\lbrack {\hat{R}\varphi} \right\rbrack_{g}^{m} = {\left( {\Sigma_{ag}^{m} + {\sum\limits_{g^{\prime} \neq g}^{G}\Sigma_{g^{\prime}g}^{m}}} \right)\varphi_{g}^{m}}} \\{\left\lbrack {\hat{S}\varphi} \right\rbrack_{g}^{m} = {\sum\limits_{g^{\prime} \neq g}^{G}{\Sigma_{{gg}^{\prime}}^{m}\varphi_{g^{\prime}}^{m}}}} \\{\left\lbrack {\hat{F}\varphi} \right\rbrack_{g}^{m} = {\chi_{g}{\sum\limits_{g^{\prime} = 1}^{G}{{\nu\Sigma}_{{fg}^{\prime}}^{m}\varphi_{g^{\prime}}^{m}}}}} \\{\left\lbrack {\hat{G}J} \right\rbrack_{gu}^{m} = {\sum\limits_{{u = x},y,z}{\frac{1}{a_{u}^{m}}{\theta_{gu}^{m}\left( {j_{gul}^{+ m} + j_{gur}^{- m}} \right)}}}}\end{matrix} \right. & (4)\end{matrix}$

An isotropic outcurrent generation operator {circumflex over (Π)} and amono-directional current throughflow operator {circumflex over (Ω)} aredefined by:

$\begin{matrix}\left\{ \begin{matrix}{\left\lbrack {\hat{\Pi}\varphi} \right\rbrack_{gu}^{m} = {c_{1{gu}}^{m}\varphi_{g}^{m}}} \\{\left\lbrack {\hat{\Omega}J} \right\rbrack_{gul}^{- m} = {{c_{2{gu}}^{m}j_{gul}^{+ m}} + {c_{3{gu}}^{m}j_{gur}^{- m}}}} \\{\left\lbrack {\hat{\Omega}J} \right\rbrack_{gur}^{+ m} = {{c_{3{gu}}^{m}j_{gul}^{+ m}} + {c_{2{gu}}^{m}j_{gur}^{- m}}}}\end{matrix} \right. & (5)\end{matrix}$

A coupling operator Ŷ couples the outcurrents to the incurrents forneighbouring cubes 10, and is defined by:

$\begin{matrix}\left\{ \begin{matrix}{j_{gul}^{+ {(m)}} = {\left\lbrack {\hat{Y}J} \right\rbrack_{gul}^{+ {(m)}} = j_{gur}^{+ {({m - 1})}}}} \\{j_{gur}^{- {(m)}} = {\left\lbrack {\hat{Y}J} \right\rbrack_{gur}^{- {(m)}} = j_{gul}^{- {({m + 1})}}}} \\{j_{gul}^{- {(m)}} = {\left\lbrack {\hat{Y}J} \right\rbrack_{gul}^{+ {(m)}} = j_{gur}^{- {({m - 1})}}}} \\{{j_{gur}^{+ {(m)}} = {\left\lbrack {\hat{Y}J} \right\rbrack_{gur}^{- {(m)}} = j_{gul}^{+ {({m + 1})}}}},}\end{matrix} \right. & (6)\end{matrix}$

This coupling operator Ŷ uses the fact that, for example, for a cube m,the u-directional left-oriented outcurrent equals the u-directionalleft-oriented incurrent for the left neighbour in direction of theCartesian axis u, and that similar equalities are verified for the otherdirections and orientations. For example, the neutron outcurrents comingfrom respective neighbours m′₁(m), m′₂(m), m′₃m), m′₄(m) into the cube mare the neutron incurrents for the cube m, as the neighbours m′₁(m),m′₂(m), m′₃m), m′₄(m) shown in FIG. 2 are respectively the neighbours(m+1) for the Cartesian axis y, (m−1) for the Cartesian axis x, (m+1)for the Cartesian axis x and (m−1) for the Cartesian axis y of Equation(6).

Using the above operators, the equations (1) and (2) can be written as:

$\begin{matrix}\left\{ \begin{matrix}{{\hat{R}\varphi} = {{\left( {{\lambda \hat{F}} + \hat{S}} \right)\varphi} + {\hat{G}j^{({in})}}}} \\{j^{({out})} = {{\hat{\Pi}\varphi} + {\hat{\Omega}j^{({in})}}}} \\{{j^{({in})} = {\hat{Y}j^{({out})}}},}\end{matrix} \right. & (7)\end{matrix}$

or as:

$\begin{matrix}{{\begin{pmatrix}{\hat{R} - \hat{S} - {\lambda \hat{F}}} & \hat{0} & \hat{G} \\{- \hat{\Pi}} & \hat{1} & {- \hat{\Omega}} \\\hat{0} & {- \hat{Y}} & \hat{1}\end{pmatrix}\begin{pmatrix}\varphi \\j^{({out})} \\j^{({in})}\end{pmatrix}} = \begin{pmatrix}0 \\0 \\0\end{pmatrix}} & (8)\end{matrix}$

which is the eigensystem to be solved.

In Equation (8), Φ is a neutron flux column vector, wherein each elementis a neutron flux Φ_(g) ^(m) for a respective cube m and for arespective energy group g (FIG. 1). Thus, the dimensions of the neutronflux column vector Φ are equal to (G×M, 1), where G is the number ofenergy groups and M is the number of cubes 10. j^((out)) is a neutronoutcurrent column vector, wherein each element is a respective neutronoutcurrent j_(gul) ^(−m), j_(gur) ^(+m), for the respective cube m,energy group g and Cartesian axis u, with u equal to x, y or z. j^((in))is a neutron incurrent column vector, wherein each element is arespective neutron incurrent j_(gul) ^(+m), j_(gur) ^(−m), for therespective cube m, energy group g and Cartesian axis u. Thus, thedimensions of the neutron outcurrent vector column j^((out)) and theneutron incurrent vector column j^((in)) are equal to (6×G×M, 1). In thefollowing of the description, the neutron outcurrent vector columnj^((out)) is also noted j_(out).

Thus, the components of an iterant (Φ, j^((out)), j^((in))) of theeigensystem defined by Equation (8) correspond to the neutron fluxesΦ_(g) ^(m), the neutron outcurrents j_(gul) ^(−m), j_(gur) ^(+m), andthe neutron incurrents j_(gul) ^(+m), i_(gur) ^(−m), to be calculatedfor each cube m and for each energy group g, and which are therespective elements of the neutron flux column vector Φ.. the neutronoutcurrent column vector j^((out)) and the neutron incurrent columnvector j^((in)).

The iterative solving procedure comprises a substep of conditioning theeigensystem into a spare eigensystem wherein the components of aniterant (j^((out)), j^((in))) of the spare eigensystem only correspondeither to the neutron incurrents j_(gul) ^(+m),j_(gur) ^(−m) coming intothe respective cube m or to the neutron outcurrents j_(gul)^(−m),j_(gur) ^(+m) coming from the respective cube m, which are therespective elements of the neutron outcurrent column vector j^((out))and the neutron incurrent column vector j^((in)).

The conditioning of the eigensystem into a spare eigensystem starts frommodifying the first part of Equation (8) written as:

({circumflex over (R)}−Ŝ−λ{circumflex over (F)})φ−Ĝj ^((in))=0,   (9)

Defining the symbolic operator {circumflex over (P)}_(λ)=({circumflexover (R)}−Ŝ−λ{circumflex over (F)})⁻¹, the neutron flux column vector Φis expressed in function of the neutron incurrent column vector j^((in))as:

φ={circumflex over (P)} _(λ) Ĝj ^((in))   (10)

By substituting this expression in the outcurrent equation part ofEquation (7), i.e. in the second part of Equation (8), a currents-onlyrelationship is obtained:

j ^((out)) =[{circumflex over (Π)}{circumflex over (P)} _(λ)Ĝ+{circumflex over (Ω)}]j ^((in))   (11)

After defining the following operator notations:

{circumflex over (Θ)}_(λ) ={circumflex over (Π)}{circumflex over (P)}_(λ)  (12)

and

{circumflex over (B)} _(λ)={circumflex over (Θ)}_(λ) Ĝ+{circumflex over(Ω)},   (13)

the currents-only relationship is written as:

j ^((out)) ={circumflex over (B)} _(λ) j ^((in))   (14)

which is a spare eigensystem.

Using the relationship j^((in))=Ŷj^((out)) between the neutronoutcurrent column vector j^((out)) and the neutron incurrent columnvector j^((in)), where Ŷ is the coupling operator, a spare eigensystemis obtained wherein the components of an iterant (j^((out))) of thespare eigensystem only correspond to the neutron outcurrents j_(gul)^(−m),j_(gur) ^(+m) coming from the respective cube m, which are theelements of the neutron outcurrent column vector j^((out)). This spareeigensystem corresponds to:

[{circumflex over (1)}−{circumflex over (B)} _(λ) Ŷ]j ^((out))=0   (15)

Since the operator {circumflex over (B)}_(λ)Ŷ will initially not beexactly equal to 1, and will become equal to 1 only as j^((out))converges to the exact solution j^((out)) _(exact) and if simultaneouslythe first neutron eigenvalue μ converges to an exact value λ_(exact), sothat Equation (15) becomes:

[{circumflex over (1)}−μ{circumflex over (B)} _(λ) Ŷ]j ^((out))=0,  (16)

with a second neutron eigenvalue μ converging towards 1, if the firstneutron eigenvalue λ converges to the exact value λ_(exact).

For numerical optimization, Equation (16) may be modified through apremultiplication with the operator Π⁻¹ and the following equation isobtained:

{circumflex over (Π)}⁻¹ j ^((out)) =μ[{circumflex over (P)} _(λ)Ĝ+{circumflex over (Π)} ⁻¹ {circumflex over (Ω)}]Ŷj ^((out))   (17)

The term {circumflex over (P)}_(λ)ĜŶj^((out)) is an isotropic term,which does not depend on the Cartesian direction u.

The expressions of the respective terms from Equation (17) for the cubeindex m, Cartesian axes u, u′ and the respective orientations s, s′along the Cartesian axes u, u′ are given by:

$\begin{matrix}\left\{ \begin{matrix}{A_{mGus} = {\sum\limits_{g \in G}\left\lbrack {{\hat{\Pi}}^{- 1}j^{({out})}} \right\rbrack_{mgus}}} \\{Q_{mG}^{G^{\prime}u^{\prime}s^{\prime}} = {\sum\limits_{g \in G}{\sum\limits_{g^{\prime} \in G^{\prime}}\left\lbrack {{\hat{P}}_{\lambda}\hat{G}\hat{Y}j^{({out})}} \right\rbrack_{mg}^{g^{\prime}u^{\prime}s^{\prime}}}}} \\{C_{mGus}^{s^{\prime}} = {\sum\limits_{g \in G}\left\lbrack {{\hat{\Pi}}^{- 1}\hat{\Omega}\hat{Y}j^{({out})}} \right\rbrack_{mgus}^{s^{\prime}}}}\end{matrix} \right. & {(18)\mspace{14mu} (19)\mspace{14mu} (20)}\end{matrix}$

G, G′ are sets of energy groups. Â and Ĉ are mono-energetic operatorswith their respective factors A_(mGus) and C_(mGus) ^(s′).Q_(mG)^(G′u′s′) are the factors of a spectral operator {circumflex over(Q)}_(λ). It should be noted that each combination {u, s}, respectively{u′, s′}, defines an interface surface of the cube m with u, u′ equal tox, y or z, and s, s′ equal to 1 or r. For example, {x, 1} defines theleft interface surface of the cube m along the Cartesian axis x.

The individual terms are given by:

$\begin{matrix}{\left\lbrack {{\hat{\Pi}}^{- 1}j^{({out})}} \right\rbrack_{mgus} = {\frac{1}{c_{1\; {gu}}^{m}}j_{mgus}^{({out})}}} & (21) \\\begin{matrix}{\left\lbrack {{\hat{P}}_{\lambda}\hat{G}\hat{Y}j^{({out})}} \right\rbrack_{mg}^{g^{\prime}u^{\prime}s^{\prime}} = {\sum\limits_{g^{\prime} = 1}^{ng}\; {{\hat{P}}_{{gg}^{\prime}}\left\lbrack {\hat{G}\hat{Y}j^{({out})}} \right\rbrack}_{mgus}^{u^{\prime}s^{\prime}}}} \\{= {\sum\limits_{g^{\prime} = 1}^{ng}\; {{\hat{P}}_{{gg}^{\prime}}\frac{1}{a_{u^{\prime}}^{m}}\theta_{g^{\prime}u^{\prime}}^{m}j_{g^{\prime}}^{m\leftarrow{\eta {({{mu}^{\prime}s^{\prime}})}}}}}}\end{matrix} & (22) \\{\left\lbrack {{\hat{\Pi}}^{- 1}\hat{\Omega}\hat{Y}j^{({out})}} \right\rbrack_{mgus}^{s^{\prime}} = {{\frac{1}{c_{1\; {gu}}^{m}}\left\lbrack {{c_{2\; {gu}}^{m}\delta_{{ss}^{\prime}}} + {c_{3\; {gu}}^{m}\left( {1 - \delta_{{ss}^{\prime}}} \right)}} \right\rbrack}j_{g}^{m\leftarrow{\eta {({mus}^{\prime})}}}}} & (23)\end{matrix}$

Equation (17) is then solved, e.g. by using a conventional iterativesolving procedure as a Gauss-Seidel procedure.

Thus, a solution is calculated for the elements j^((out)) _(mgus) of theneutron outcurrent column vector j^((out)). This enables to determinethe solution of the neutron incurrent column vector j^((in)) accordingto Equation (14), and finally to determine the solution of the neutronflux column vector Φ according to Equation (10).

The calculated neutron flux column vector Φ obtained through themodelizing method can be used to control an existing nuclear reactorcore, e.g. for managing the nuclear fuel, or be used for building a newreactor core.

The modelizing method may be implemented on parallel processors or on asingle processor.

This above-disclosed modelizing method has proved to lead to betterconvergence accuracy, a better computational robustness and a bettercomputational efficiency. This is connected to the solving of a sparseeigensystem wherein the only components of the iterant correspond toneutron outcurrents j^((out)), and do not depend on neutron fluxes.

Further, the convergence accuracy of the modelizing method according tothe invention may be improved up to 1E-12, whereas the convergenceaccuracy obtained with a classic modelizing method is limited to 1E-6.

In other embodiments, the components of the iterant of the spareeigensystem to be solved correspond only to neutron incurrents j^((in)).

In order to further improve robustness and computational efficiency,according to a second aspect of the invention, the eigensystem is firstconditioned in a restricted eigensystem corresponding to the eigensystemfor a selection of some neutron energy groups. This selection is alsocalled spectral restriction to some energy groups. The restrictedeigensystem may then be solved according to the first aspect of theinvention. Finally, the solution of the restricted eigensystem is usedto solve the eigensystem.

The number NGC of the selected energy groups is smaller than the totalnumber ng of energy groups. NGC may be the number of coarsened spectralbands which are collections of fine energy groups merged into a smallernumber of coarse energy groups.

ng is, for example, equal to 8, and NGC may for example be equal to 4,3, 2 or 1.

According to this second aspect, the Equation (17) for the drivingfactors d_(Gmus) ^((out)), with its respective terms given by Equations(18) to (20), is expressed as:

$\begin{matrix}{\underset{\underset{{}_{}^{}{}_{}^{}}{}}{A_{Gmus}^{({out})}d_{Gmus}^{({out})}} = {\mu\left\lbrack {\underset{\underset{{\,^{``}{isotropic}}\mspace{14mu} {production}^{''}}{}}{\sum\limits_{G^{\prime} = 1}^{NGC}\; {\sum\limits_{u^{\prime} = 1}^{3}\; {\sum\limits_{{s^{\prime} = 1},r}\; {Q_{mG}^{G^{\prime}u^{\prime}s^{\prime}}d_{G^{\prime}{\eta {({{mu}^{\prime}s^{\prime}})}}u^{\prime}s^{\prime}}^{({out})}}}}} + \underset{\underset{{}_{}^{}{}_{}^{}}{}}{\sum\limits_{{s^{\prime} = 1},r}\; {C_{mGus}^{s^{\prime}}d_{{Gmus}^{\prime}}^{({out})}}}} \right\rbrack}} & (24)\end{matrix}$

where G, G′ are coarse neutron energy group indexes, and

in which an isotropic term

$\sum\limits_{G^{\prime} = 1}^{NGC}\; {\sum\limits_{u^{\prime} = 1}^{3}\; {\sum\limits_{{s^{\prime} = 1},r}\; {Q_{mG}^{G^{\prime}u^{\prime}s^{\prime}}d_{G^{\prime}{\eta {({{mu}^{\prime}s^{\prime}})}}u^{\prime}s^{\prime}}^{({out})}}}}$

is independent of the outgoing current direction specified by theCartesian axis u and orientation s, and a throughflow term

${\sum\limits_{{s^{\prime} = 1},r}\; {C_{mGus}^{s^{\prime}}d_{{Gmus}^{\prime}}^{({out})}}},$

also called anisotropic term, is defined within the same coarse neutronenergy group G along the same Cartesian axis u.

The decomposition of the outcurrents leaving the cube m into isotropicand anisotropic terms is illustrated on FIG. 3.

The isotropic term is the same in all directions and the anisotropicterm varies between interface surfaces of the cube m.

Solving at first the eigensystem for a spectral restriction to someenergy groups, moreover with a decomposition of the outcurrents in whichthe isotropic term is computed easily, leads to a better computationalrobustness and a better computational efficiency. This is connected tothe reduction of eigensystem dimensions that results from the spectralrestriction according to this second aspect of the invention.

In order to further improve robustness and computational efficiency,according to a third aspect of the invention, the iterative solvingprocedure for a plurality of energy groups may be in form of amulti-level V-cycle 20, as shown on FIG. 4, comprising a top level 21, afirst intermediate level 22, four second intermediate levels 24, 26, 28,30, and three bottom levels 32, 34, 36.

The top level 21 of the V.-cycle comprises the iteration for theeigensystem corresponding to the steady-state diffusion equation (1), orNEM equation, for the plurality of neutron energy groups, and theconditioning of the eigensystem into a restricted eigensystem for aspectral restriction to some energy groups according to the secondaspect of the invention. The restricted eigensystem resulting from thetop level 21 is then fed into a first intermediate level 22 just underthe top level 21.

The first intermediate level 22 comprises the conditioning of therestricted eigensystem into a spare eigensystem according to the firstaspect of the invention wherein the components of an iterant (j^((out)),j^((in))) of the spare eigensystem only correspond either to neutronincurrents j_(gul) ^(+m), j_(gur) ^(−m) coming into the cubes m or toneutron outcurrents j_(gul) ^(−m), j_(gur) ^(+m), also noted j^((out))_(mgus), coming from the cubes m. The factors of operators Â,{circumflex over (Q)}_(λ) and Ĉ given by Equations (18) to (20) arecomputed explicitly for this first intermediate level 22.

Each second intermediate levels 24, 26, 28, 30 comprises theconditioning of a former spare eigensystem for a former selection ofneutron energy groups into a latter spare eigensystem for a latterselection of neutron energy groups, the number of neutron energy groupsin said latter selection being smaller to the number of neutron energygroups in said former selection. The former spare eigensystem of thesecond intermediate level 24 subsequent to the first intermediate level22 in the downward orientation of the multi-level V-cycle 20 is thespare eigensystem resulting from the first intermediate level 22. Theformer spare eigensystem of each second intermediate level 26, 28, 30which is not subsequent to the first intermediate level 22 correspondsto the latter spare eigensystem resulting from the precedent secondintermediate level 24, 26, 28. In other words, the number of neutronenergy groups decreases from a second intermediate level to the nextsecond intermediate level in the downward orientation. At the lastsecond intermediate level 30, the number of neutron energy groups may beequal to 1. It should be noted that this restriction of the number ofneutron energy groups for the second intermediate levels 24, 26, 28, 30may be done according to a spectral condensation algebra, which will bedescribed later, and not according to the second aspect of theinvention.

Bottom levels 32, 34, 36 comprise, in the downward orientation of theV-cycle, the solving, according to state of the art procedures, of thelast spare eigensystem for a single energy group determined at the lastsecond intermediate level 30. Bottom levels 34, 36 correspond to aspatial restriction of the eigensystem resulting from the precedentbottom level 32, 34 in respective coarsened grids with cubes beinglarger than the cubes of the precedent bottom level 32, 34.

At the bottom level 36 of the V-cycle, a solution of the eigensystem forthe single energy group is computed.

This solution is then reintroduced in the upper levels of the V-cycle inthe upward orientation so that solutions are computed for the respectivespare eigensystems.

At the top level 21 of the V-cycle in the upward orientation, thesolution of the NEM equation for the plurality of energy groups iscomputed.

In the illustrated embodiment of FIG. 4, the number ng of neutron energygroups for the NEM equation at the top level 21 may be equal to 8. Forthe intermediate levels 22, 24, 26, 28, 30, the number of neutron energygroups in the respective selections may be respectively equal to four,three, two and one. The three bottom levels 32, 34, 36 of the V-cycleshown on FIG. 4 comprise the solving of said last spare eigensystem forsaid one energy group according to the CMR procedure.

The four intermediate levels 22, 24, 26, 28 corresponding to the solvingof respective spare eigensystems are also designated by the respectivereferences SR4, SR3, SR2 and SR1 (FIG. 5).

The factors of the mono-energetic operator Â for the level SR4 arecomputed explicitly. Subsequently, the factors of the mono-energeticoperator Â for the level SR3 are derived from the ones for the levelSR4, as shown on FIG. 5, where each box 37 represents a neutron energygroup, through:

A_(m;SR3) ^(1u′s′)=A_(m;SR4) ^(1u′s′)

A_(m;SR3) ^(2u′s′)=A_(m;SR4) ^(2u′s′)

A _(m;SR3) ^(3u′s′) =A _(m;SR4) ^(3u′s′) +A _(m;SR4) ^(4u′s′)  (25)

Equation (25) is illustrated on FIG. 5 by the arrows 38 between levels22 and 24.

Subsequently, the factors of the mono-energetic operator Â for the levelSR2 are derived from the ones for the level SR3 through:

A _(m;SR2) ^(1u′s′) =A _(m;SR3) ^(1u′s′) +A _(m;SR3) ^(2u′s′)

A_(m;SR2) ^(2u′s′)=A_(m;SR3) ^(3u′s′)  (26)

Equation (26) is illustrated on FIG. 5 by the arrows 40 between levels24 and 26.

Finally, the factors of the mono-energetic operator Â for the level SRIare derived from the ones for the level SR2 through:

A _(m;SR1) ^(u′s′) =A _(m;SR2) ^(1u′s′) +A _(m;SR2) ^(2u′s′)  (27)

Equation (27) is illustrated on FIG. 5 by the arrows 42 between levels26 and 28.

Equations (25) to (27) define a spectral condensation algebra for themono-energetic operator Â corresponding to a spectral grid partitioningstrategy shown in FIG. 5. It should be noted that another spectralcondensation algebra for the mono-energetic operator Â with differentequations may be defined with another spectral grid partitioningstrategy, i.e. with a different arrangement of the boxes 37.

The factors of the mono-energetic operator Ĉ are derived in a similarmanner as for the factors of the mono-energetic operator Â, from thelevel SR4 to the level SR1 in the downward orientation.

The factors of the spectral operator {circumflex over (Q)}_(λ) for thelevel SR4 are computed explicitly. Subsequently, the factors of thespectral operator {circumflex over (Q)}_(A) for the level SR3 arederived from the ones for the level SR4 through:

Q_(m1;SR3) ^(1u′s′)=Q_(m1;SR4) ^(1u′s′)

Q_(m2;SR3) ^(2u′s′)=Q_(m2;SR4) ^(2u′s′)

Q_(m1;SR3) ^(2u′s′)=Q_(m1;SR4) ^(2u′s′)

Q_(m2;SR3) ^(1u′s′)=Q_(m2;SR4) ^(1u′s′)

Q _(m2;SR3) ^(3u′s′) =Q _(m2;SR4) ^(3u′s′) +Q _(m2;SR4) ^(4u′s′)

Q _(m3;SR3) ^(2u′s′) =Q _(m3;SR4) ^(2u′s′) +Q _(m4;SR4) ^(2u′s′)

Q _(m1;SR3) ^(3u′s′) =Q _(m1,SR4) ^(3u′s′) +Q _(m1;SR4) ^(4u′s′)

Q _(m3;SR3) ^(1u′s′) =Q _(m3;SR4) ^(1u′s′) +Q _(m4;SR4) ^(1u′s′)

Q _(m3;SR3) ^(3u′s′) =Q _(m3;SR4) ^(3u′s′) +Q _(m3;SR4) ^(4u′s′) +Q_(m4;SR4) ^(3u′s′) +Q _(m4;SR4) ^(4u′s′)  (28)

Subsequently, the factors of the spectral operator {circumflex over(Q)}_(λ) for the level SR2 are derived from the ones for the level SR3through:

Q _(m2;SR2) ^(1u′s′) =Q _(m3;SR3) ^(1u′s′) +Q _(m3;SR3) ^(2u′s′)

Q _(m1;SR2) ^(2u′s′) =Q _(m1;SR3) ^(3u′s′) +Q _(m2;SR3) ^(3u′s′)

Q_(m2;SR2) ^(2u′s′)=Q_(m2;SR3) ^(2u′s′)

Q _(m1;SR2) ^(1u′s′) =Q _(m1;SR3) ^(1u′s′) +Q _(m1;SR3) ^(2u′s′) +Q_(m2;SR3) ^(1u′s′) +Q _(m2;SR3) ^(2u′s′)  (29)

Finally, the factors of the spectral operator {circumflex over (Q)}_(λ)for the level SR1 are derived from the ones for the level SR2 through:

Q _(m;SR1) ^(u′s′) =Q _(m1;SR2) ^(1u′s′) +Q _(m1;SR2) ^(2u′s′) +Q_(m2;SR2) ^(1u′s′) +Q _(m2;SR2) ^(2u′s′)  (30)

Equations (28) to (30) define a spectral condensation algebra for thespectral operator {circumflex over (Q)}_(λ) corresponding to a spectralgrid partitioning strategy not shown and similar to the one shown inFIG. 5 for the mono-energetic operators Â and Ĉ. It should be noted thatanother spectral condensation algebra for the spectral operator{circumflex over (Q)}_(λ) may also be define with another spectral gridpartitioning strategy.

Thus, a solution is calculated for the elements j_((out)) _(mgus) of theneutron outcurrent column vector j^((out)) at the top level 21 of theV-cycle in the upward orientation. This enables to determine thesolution of the neutron incurrent column vector j^((in)) according toEquation (14), and finally to determine the solution of the neutron fluxcolumn vector Φ according to Equation (10).

The nuclear reactor core is then built or operated on the basis of thecalculated neutron flux column vector Φ.

Using spectral condensation algebras, wherein no complex arithmeticoperation is involved, for the operators Â, Ĉ and {circumflex over(Q)}_(λ) at the respective levels SR3, SR2 and SR1 allows the operatorsÂ, Ĉ and {circumflex over (Q)}_(λ) to be computed very cheaply, thusleading to a better computational robustness and a better computationalefficiency of the modelizing method.

In order to further improve robustness and computational efficiency,according to a fourth aspect of the invention, the cubes 10 are splitinto a first category and a second category.

In the following description and as shown on FIG. 6, the cubes 10R ofthe first category will be called the red cubes and the cubes 10B of thesecond category will be called the black cubes but no specificrestrictive meaning should be associated with the words “black” and“red”. Each red cube 10R has only black cubes 10B as direct neighbours.Thus, most of the red cubes 10R have six direct black neighbours 10B. Itshould be understood by “direct” neighbours, the cubes sharing a commonsurface with the considered cube.

Consequently, and as illustrated by FIG. 6, the grid 12 has a visualanalogy with respect to the dark and light regions of a checkerboard ina two-dimensional representation.

Then, the cubes are numbered, starting for example by the red cubes 10Rand ending by the black cubes 10B.

In the following description, such a split of the cubes in twocategories and the numbering of one category after the other will bereferred to as red-black ordering.

An advantage of the red-black ordering of the cubes in comparison withthe state of the art lexicographical ordering is that, for a red cube1R, all its direct neighbours will be black, and vice versa.

The Equation (24), which the computer has to solve in order to calculateneutron flux within the core, can be written in the matrix-vector form:

Âd=μ[{circumflex over (Q)} _(λ) +Ĉ]d   (31)

Using a variable parameter γ, which value is typically comprised between0.9 and 0.95, so that the product γμ forms a shift, an iteration ofEquation (31) is written with the following shift-inverted implicitform:

[Â−γμ ^((n−1))({circumflex over (Q)} _(λ+Ĉ))] d ^((n)) =s ^((n−1))  (32)

where s^((n−1)) is a source given by:

s ^((n−1))=(1−γ)μ^((n−1))({circumflex over (Q)} _(λ+Ĉ)) d ^((n−1))  (33)

The spectral operator {circumflex over (Q)}_(λ) and the mono-energeticoperator Ĉ are sparse, coupling red cubes 10R only to direct blackneighbours 10B and vice versa. Thus, the red-black ordering enables thefollowing convenient relationship between the red and the black parts ofthe equation:

Âd _(red)−γμ^((n−1))({circumflex over (Q)} _(λ) +Ĉ) d _(black) =s _(red)

Âd _(black)−γμ^((n−1))({circumflex over (Q)} _(λ) +Ĉ) d _(red) =s_(black)   (34)

Equation (34) is transformed into the following equivalent equation forthe red solution part only:

$\begin{matrix}{{{\overset{\sim}{\Theta}}_{red}{\underset{\_}{d}}_{red}} = {{\overset{\sim}{\underset{\_}{s}}}_{red}\mspace{14mu} {with}\mspace{14mu} \left\{ \begin{matrix}{{\overset{\sim}{\Theta}}_{red} = {\hat{1} - {\left( {\gamma\mu}^{({n - 1})} \right)^{2}\left\lbrack {{\hat{A}}^{- 1}\left( {{\hat{Q}}_{\lambda} + \hat{C}} \right)} \right\rbrack}^{2}}} \\{{{\overset{\sim}{\underset{\_}{s}}}_{red} = {{\hat{A}}^{- 1}\begin{bmatrix}{{\underset{\_}{s}}_{red} + {{\gamma\mu}^{({n - 1})}{\hat{A}}^{- 1}}} \\{\left( {{\hat{Q}}_{\lambda} + \hat{C}} \right){\hat{A}}^{- 1}{\underset{\_}{s}}_{black}}\end{bmatrix}}},}\end{matrix} \right.}} & (35)\end{matrix}$

From a calculation point of view, the use of a red-black ordering meansthat, during an iterative solving procedure if, in the first half of aniteration, the red components d _(red) of the iterant d are updated,then, during the second part of the iteration, all the black componentsd _(black) will be updated on the basis of the red components of theirred neighbours d _(red). In other words, the value for each black cube10B will be calculated on the basis of the values for all its direct redneighbours 10R.

Using the red-black ordering improves the computation of the operatorsÂ, Ĉ and {circumflex over (Q)}_(λ) described in the first, second orthird aspects of the invention. Thus, the red-black ordering may be usedin complement of the first, second or third aspects of the invention inorder to improve the computational efficiency.

Such a red-black ordering proves to be especially useful when used witha particular solving procedure which constitutes a fifth aspect of theinvention.

This procedure is a particular highly robust stabilized Bi-ConjugateGradient Stabilized (Bi-CGStab) procedure. An adequate introductorydescription of this procedure can be found in the following references:

-Y. Saad, “Iterative Methods for Sparse Linear Systems”, second edition,Society for Industrial and Applied Mathematics (SIAM) (2003);

H. A. van der Vorst, “Bi-CGSTAB: a Fast and Smoothly Converging Variantof Bi-CG for the solution of nonsymmetric linear systems”, SIAM J. Sci.Stat. Comput. 13(2), pp. 631-644 (1992),

The Bi-CGStab procedure is derived from a minimization principle for afunctional of d, with given Θ and s, for which stationarity applies withregard to small variations δd around the specific optimum d whichsatisfies, exactly, the linear system given by Equation (35), and forwhich the functional assumes a minimum value.

Thus, it is possible to define a solving procedure driven by theminimization of a functional rather than by more direct considerationson how to solve Equation (35) efficiently. The Bi-CGStab procedurefollows from such a minimization principle, where the successive changesin the iterant are organized such that each change in the functional(which is like a Galerkin-weighted integrated residual) is orthogonalwith respect to all of the preceding changes.

The particular Bi-CGStab procedure according to the fifth aspect of theinvention is given below with solution vector d, solution residual r(with r=s−Θd) and auxiliary vector r*, s and p, and with initialsolution estimate d₀:

$\begin{matrix}\left\{ \begin{matrix}1. & {{r_{0}:={\overset{\_}{s} - {\overset{\bullet}{\Theta}d_{0}}}},{r^{*} = r_{0}}} \\2. & {{p_{0}:=r_{0}},} \\3. & {{{{do}\; i} = 0},1,\ldots \mspace{14mu},N} \\4. & {{\alpha_{i}:=\frac{\left( {r^{*},r_{i}} \right)}{\left( {r^{*},{\overset{\bullet}{\Theta}p_{i}}} \right)}},} \\5. & {{s_{i}:={r_{i} - {\alpha_{i}\overset{\bullet}{\Theta}p_{i}}}},} \\6. & {{\omega_{i}:=\frac{\left( {{\overset{\bullet}{\Theta}\; s_{i}},s_{i}} \right)}{\left( {{\overset{\bullet}{\Theta}\; s_{i}},{\overset{\bullet}{\Theta}s_{i}}} \right)}},} \\7. & {{d_{i + 1}:={d_{i} + {\alpha_{i}p_{i}} + {\omega_{i}s_{i}}}},} \\8. & {{r_{i + 1}:={s_{i} - {\omega_{i}\overset{\bullet}{\Theta}s_{i}}}},} \\9. & {{\beta_{i}:=\frac{\left( {r^{*},r_{i + 1}} \right)\alpha_{j}}{\left( {r^{*},r_{i}} \right)\omega_{j}}},} \\10. & {p_{i + 1}:={r_{i + 1} + {\beta_{i + 1}\left( {p_{i} - {\omega_{i}\overset{\bullet}{\Theta}p_{i}}} \right)}}} \\11. & {{end}\mspace{14mu} {do}}\end{matrix} \right. & (36)\end{matrix}$

The above Bi-CGStab sequence is truncated after N steps, with usuallyN<3.

For the Bi-CGStab procedure, the preconditioning is based on theshift-inverted implicit form of Equation (32), with typically0.9<γ<0.95.

With this choice for the shift γμ, the operator s is guaranteed toremain nonsingular since, during the solution process, γμ will convergedown to γ times the smallest possible eigenvalue of the system.

Solving Equation (32) is numerically equivalent to determining the fluxdistribution in a slightly subcritical system with a neutron source,which is a perfect scenario for deploying advanced preconditioned Krylovprocedure.

Using the red-black ordering, Equation (32) is transformed into Equation(35) as above explained. The preconditioned system given by Equation(35) constitutes the sixth aspect of the invention.

Once Equation (35) has been solved through the Bi-CGStab procedure theblack solution part is to be computed from the black part of Equation(34).

Use of this particular way of preconditioning, which means that thedirect neutron interaction between neighboring (red vs. black) nodes, asprojected onto the grid, is preincluded in the system to be solved,making the unity operator in Equation (35) more dominant since

∥(Â ⁻¹({circumflex over (Q)} _(λ) +Ĉ))² ∥<∥Â ⁻¹({circumflex over (Q)}_(λ) +Ĉ)∥,   (37)

given that

∥Â ⁻¹({circumflex over (Q)} _(λ) +Ĉ)∥<1   (38)

This way of preconditioning manages to pre-include, at low computationalcost, crucial information (or the major part of that information) withregard to the physical interactions between nodes that determine thespatial couplings and hence the solution of the equation.

In another embodiment, the iterative solving procedure is a Gauss-Seidelprocedure, or an iterative procedure in the Krylov subspace, such as ageneralized minimal residual method, also abbreviated GMRES, or a plainJacobi method.

According to a seventh aspect of the invention, a variational principlefor improved eigenvalue computation will be described in the followingdescription.

Equation (16) corresponding to the spare eigensystem, above describedfor the first, second or third aspect of the invention, is written inthe following form, with the second neutron eigenvalue μ and the neutronoutcurrent vector column j_(out):

j _(out) =μ{circumflex over (Π)}{circumflex over (P)} _(λ) ĜŶj _(out)+{circumflex over (Ω)}Ŷj _(out),

j _(out) =μc ₁ φ+{circumflex over (Ω)}Ŷj _(out)   (39)

where c₁ is another notation for the operator {circumflex over (Π)}depending on the cube properties.

In what follows, we will assume that the first neutron eigenvalue λfulfills the role of a scalar control parameter. Using the property thatthe second neutron eigenvalue μ approaches unity if the total solutionconverges and thus stabilizes, the variational principle, orperturbation approach, comprises the varying the first neutroneigenvalue λ to drive the second neutron eigenvalue μ towards 1 and toaccelerate the convergence of the spare eigensystem solution. Thisvariation of the first neutron eigenvalue λ, may be done prior to everyiteration of the spare eigensystem, or prior to every second, third,fourth or fifth iteration of the spare eigensystem.

The variation of the first neutron eigenvalue λ is an application of avariation δλ to λ such that λ is driven to λ+δλ, which impacts thesecond neutron eigenvalue μ through a perturbed interface currentequation:

$\begin{matrix}\begin{matrix}{{j_{out} + {\delta \; j_{out}}} = {\left\lbrack {{\left( {\mu + {\delta\mu}} \right){\hat{\Pi}\left( {{\hat{P}}_{\lambda} + {{\delta\lambda}\frac{\partial{\hat{P}}_{\lambda}}{\partial\lambda}}} \right)}\hat{G}} + \hat{\Omega}} \right\rbrack {\hat{Y}\left( {j_{out} + {\delta \; j_{out}}} \right)}}} \\{= {{\left( {\mu + {\delta\mu}} \right){c_{1}\left( {\varphi + {{\delta\lambda}\frac{\partial\varphi}{\partial\lambda}}} \right)}} + {\hat{\Omega}{\hat{Y}\left( {j_{out} + {\delta \; j_{out}}} \right)}}}}\end{matrix} & (40)\end{matrix}$

with the partial operator derivative:

$\begin{matrix}{{\hat{\Pi}\frac{\partial{\hat{P}}_{\lambda}}{\partial\lambda}\hat{G}\hat{Y}j_{out}} \equiv {c_{1}\frac{\partial\varphi}{\partial\lambda}}} & (41)\end{matrix}$

As perturbations δj_(out) are expected to be increasingly small duringthe iterative procedure, they are conveniently ignored in Equation (40),and the integration with an arbitrary adjoined weighting functionω^(†)yields the following expression:

$\begin{matrix}{{{\delta\mu}\lbrack{\delta\lambda}\rbrack} \cong {\frac{\langle{\omega^{\dagger}{c_{1}\varphi}}\rangle}{\langle{\omega^{\dagger}{c_{1}\frac{\partial\varphi}{\partial\lambda}}}\rangle}{\delta\lambda}}} & (42)\end{matrix}$

As the variation of the first neutron eigenvalue λ is such that thesecond neutron eigenvalue μ is driven towards 1, i.e. that μ+δμ, isdriven towards 1, the approximation δμ≡1−μ is imposed. With thisapproximation and the Equation (42), the following equation is obtained:

$\begin{matrix}{{\delta\lambda} \cong {\frac{\langle{\omega^{\dagger}{c_{1}\varphi}}\rangle}{\langle{\omega^{\dagger}{c_{1}\frac{\partial\varphi}{\partial\lambda}}}\rangle}\left( {1 - \mu} \right)}} & (43)\end{matrix}$

Using the equivalence μc₁φ=j_(out)−{circumflex over (Ω)}Ŷj_(out), δλ isgiven by:

$\begin{matrix}{{{\delta\lambda} \cong \frac{\langle{\omega^{\dagger}{j_{out} - {c_{1}\varphi} + {\hat{\Omega}\hat{Y}j_{out}}}}\rangle}{\langle{\omega^{\dagger}{c_{1}\frac{\partial\varphi}{\partial\lambda}}}\rangle}} = \frac{\langle{\omega^{\dagger}r_{out}}\rangle}{\langle{\omega^{\dagger}{c_{1}\frac{\partial\varphi}{\partial\lambda}}}\rangle}} & (44)\end{matrix}$

where r_(out) represents the residual of the interface outcurrentbalance equation, i.e. Equation (39), with the imposed target value forthe second neutron eigenvalue μ equal to 1. The residual r_(out) needsto be computed for each NEM iteration step in any case.

The first neutron eigenvalue λ, considered in this aspect as a scalarcontrol parameter, is then updated such that:

$\begin{matrix}{{\lambda^{({new})} = {\lambda^{({prev})} + \frac{\langle{\omega^{\dagger}r_{out}^{({prev})}}\rangle}{\langle{\omega^{\dagger}{c_{1}\frac{\partial\varphi}{\partial\lambda}}}\rangle}}},} & (45)\end{matrix}$

where λ^((prev)) and λ^((new)) are respectively the previous and theupdated value of the first neutron eigenvalue λ.

For determining the partial operator derivative

$\frac{\partial\varphi}{\partial\lambda}$

and using Equation (7) for example, the neutron flux Φ is written:

φ=({circumflex over (R)}−Ŝ−λ{circumflex over (F)})⁻¹ ĜŶj _(out)   (46)

from which it follows that:

$\begin{matrix}{\frac{\partial\varphi}{\partial\lambda} = {\left\lbrack {\frac{\partial}{\partial\lambda}\left( {\hat{R} - \hat{S} - {\lambda \hat{F}}} \right)^{- 1}} \right\rbrack \hat{G}\hat{Y}j_{out}}} & (47)\end{matrix}$

An approximate derivative

$\frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}$

is then computed by neglecting the upper-diagonal part of the matrixrelative to the inscatter operator Ŝ such that Ŝ=Ŝ_(LD), i.e. byneglecting the upscattering:

$\begin{matrix}{\frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} = {\left\lbrack {\frac{\partial}{\partial\lambda}\left( {\hat{R} - {\hat{S}}_{LD} - {\lambda \hat{F}}} \right)^{- 1}} \right\rbrack \hat{G}\hat{Y}j_{out}}} & (48)\end{matrix}$

It should be noted that in many computational cases this is not even anapproximation, but numerically exact if no upscattering is modeled.

It is assumed that:

$\begin{matrix}{\frac{\partial\varphi}{\partial\lambda} \cong \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}} & (49)\end{matrix}$

which will be sufficient for a successful application of the variationalprinciple.

For computing

$\frac{\partial\overset{\_}{\varphi}}{\partial\lambda},$

it is first written:

{circumflex over (R)}−Ŝ _(LD) −λ{circumflex over (F)}=({circumflex over(R)}−Ŝ _(LD))({circumflex over (1)}−λ({circumflex over (R)}−Ŝ _(LD))⁻¹{circumflex over (F)})   (50)

from which it follows that:

({circumflex over (R)}−Ŝ _(LD) −λ{circumflex over (F)})⁻¹=({circumflexover (1)}−λ({circumflex over (R)}−Ŝ _(LD))⁻¹ {circumflex over(F)})⁻¹({circumflex over (R)}−Ŝ _(LD))⁻¹   (₅₁)

Then defining an operator notation {circumflex over (B)} with:

{circumflex over (B)}=({circumflex over (R)}−Ŝ _(LD))⁻¹ {circumflex over(F)}  (52)

and in order to implicitly compute

${\frac{\partial}{\partial\lambda}\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack}^{- 1},$

the Taylor formula is applied to {circumflex over (1)}−λ{circumflex over(B)}:

$\begin{matrix}{\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack^{- 1} = {{\hat{1} + {\lambda \; \hat{B}} + {\lambda^{2}{\hat{B}}^{2}} + {\lambda^{3}{\hat{B}}^{3}} + \ldots} = {\sum\limits_{n = 0}^{\infty}\left\lbrack {\lambda \; \hat{B}} \right\rbrack^{n}}}} & (53)\end{matrix}$

The Taylor expansion is rearranged as:

{circumflex over (1)}+[{circumflex over (1)}+λ{circumflex over (B)}+λ ²{circumflex over (B)} ² +λ ³ {circumflex over (B)} ³ + . . .]λ{circumflex over (B)}=[{circumflex over (1)}−λ{circumflex over (B)}]⁻¹ λ{circumflex over (B)}  (54)

By application of the chain rule to Equation (54), the derivative

${\frac{\partial}{\partial\lambda}\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack}^{- 1}$

is given by:

$\begin{matrix}{{\frac{\partial}{\partial\lambda}\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack}^{- 1} = {{\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack^{- 1}\hat{B}} + {\left\lbrack {\frac{\partial}{\partial\lambda}\left\lbrack {1 - {\lambda \; \hat{B}}} \right\rbrack}^{- 1} \right\rbrack \lambda \; \hat{B}}}} & (55)\end{matrix}$

which yields:

$\begin{matrix}{{\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack {\frac{\partial}{\partial\lambda}\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack}^{- 1}} = {\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack^{- 1}\hat{B}}} & (56)\end{matrix}$

The final expression that allows an efficient iterative scheme fordetermining the derivative

${\frac{\partial}{\partial\lambda}\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack}^{- 1}$

is obtained by premultiplication of Equation (56) with {circumflex over(1)}−λ{circumflex over (B)}, which gives:

$\begin{matrix}{{\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack^{2}{\frac{\partial}{\partial\lambda}\left\lbrack {\hat{1} - {\lambda \; \hat{B}}} \right\rbrack}^{- 1}} = \hat{B}} & (57)\end{matrix}$

By writing [{circumflex over (1)}−λ{circumflex over (B)}]²={circumflexover (1)}−2λ{circumflex over (B)}+λ²{circumflex over (B)}², computationof

$\frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}$

is obtained from:

$\begin{matrix}{{\left\lbrack {\hat{1} - {2\lambda \; \hat{B}} + {\lambda^{2}{\hat{B}}^{2}}} \right\rbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}} = s} & (58)\end{matrix}$

with the source term s defined by:

s={circumflex over (B)}({circumflex over (R)}−Ŝ _(LD))⁻¹ ĜŶj_(out)=({circumflex over (R)}−Ŝ _(LD))⁻¹ {circumflex over(F)}({circumflex over (R)}−Ŝ _(LD))⁻¹ ĜŶj _(out)   (59)

In the case of two energy groups, Equation (58) is solved directly andanalytically.

In the case of more than two energy groups, Equation (58) is solvediteratively, without exactness required in early iteration phases,through application of:

$\begin{matrix}{{\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} \right\rbrack^{({new})} = {s + {2\lambda \; {\hat{B}\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} \right\rbrack}^{({old})}} - {\lambda^{2}{B^{2}\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} \right\rbrack}^{({old})}}}},} & (60)\end{matrix}$

where

$\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} \right\rbrack^{({old})}\mspace{14mu} {{and}\mspace{14mu}\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} \right\rbrack}^{({new})}$

are respectively the previous and updated values of the approximatederivative

$\frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}.$

Early means about the first third up to the first half of the totallyneeded number of iteration steps.

It should be noted that the iterations from the computation of

$\frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}$

can be abort rather early, since this variational approach needs only anapproximate value of

$\frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}$

to be useful.

Thus, this variational principle according to this seventh aspect of theinvention is very effective, because the numerator

ω^(†)|r_(out) ^((prev))

in Equation (45) becomes smaller and smaller upon convergence, whereasthe denominator

$\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\lambda}} \right.\rangle$

will converge rather quickly to an adequate value. In this seventhaspect, the numerator defines the driving principle since it convergestowards zero.

Varying the first neutron eigenvalue λ according to Equation (45) drivesthe second neutron eigenvalue μ towards 1 and accelerates theconvergence of the spare eigensystem solution. This variation of thefirst neutron eigenvalue λ may be done prior to every iteration of thespare eigensystem, or prior to every second, third, fourth or fifthiteration of the spare eigensystem.

Thus, the variational principle of the first neutron eigenvalue λ may beused in complement of the first, second or third aspects of theinvention in order to improve the computational efficiency.

All the above mentioned products of this seventh aspect can be reducedto the black part of the grid 12 without losing numerical efficiency inthe method, so that the red-black ordering above described in the fourthaspect of the invention may be used in complement of this seventhaspect.

According to an eighth aspect of the invention, the above describedvariational principle for the eigenvalue computation can be generalizedfor the computation of a control parameter being noted Θ, that typicallymodulates the removal operator {circumflex over (R)} given by Equation(4), the modulated removal operator being noted {circumflex over(R)}_(Θ).

In the example of a modulation in all positions, the control parameter Θis the boron concentration. In the example of a modulation in selectedpositions, the control parameter Θ is a parameter for controlling a rodinsertion depth for a group of selected control rods.

It is assumed that {circumflex over (R)}={circumflex over (R)}_(Θ)

Equation (16) is then written in the following form, with the secondneutron eigenvalue μ:

j _(out) =μ{circumflex over (Π)}{circumflex over (P)} _(Θ) ĜŶj _(out)+{circumflex over (Ω)}Ŷj _(out)

j _(out) =μc ₁ φ+{circumflex over (Ω)}Ŷj _(out)   (61)

with a symbolic spectral operator {circumflex over (P)}_(Θ) defined by

{circumflex over (P)} _(Θ)=({circumflex over (R)} _(Θ) −Ŝ−{circumflexover (F)})⁻¹   (62)

Using the property that the second neutron eigenvalue μ approaches unityif the total solution, including the correct value of the controlparameter Θ, converges and thus stabilizes, the variational principle,or perturbation approach, comprises varying the control parameter Θ todrive the second neutron eigenvalue μ towards l.

The variation of the control parameter Θ is an application of avariation δΘ to Θ such that Θ is driven to Θ+δΘ, which impacts thesecond neutron eigenvalue μ through a perturbed interface currentequation:

$\begin{matrix}{{{j_{out} + {\delta \; j_{out}}} = {\left\lbrack {{\left( {\mu + {\delta \; \mu}} \right){\hat{\Pi}\left( {{\hat{P}}_{\vartheta} + {\delta \; \vartheta \frac{\partial{\hat{P}}_{\vartheta}}{\partial\vartheta}}} \right)}\hat{G}} + \hat{\Omega}} \right\rbrack {\hat{Y}\left( {j_{out} + {\delta \; j_{out}}} \right)}}}{{j_{out} + {\delta \; j_{out}}} = {{\left( {\mu + {\delta\mu}} \right){c_{1}\left( {\varphi + {\delta \; \vartheta \frac{\partial\varphi}{\partial\vartheta}}} \right)}} + {\hat{\Omega}\; {\hat{Y}\left( {j_{out} + {\delta \; j_{out}}} \right)}}}}} & (63)\end{matrix}$

with the partial operator derivative:

$\begin{matrix}{{\hat{\Pi}\frac{\partial{\hat{P}}_{\vartheta}}{\partial\vartheta}\hat{G}\hat{Y}j_{out}} \equiv {c_{1}\frac{\partial\varphi}{\partial\vartheta}}} & (64)\end{matrix}$

As perturbations δj_(out) are expected to be increasingly small duringthe iterative procedure, they are conveniently ignored in Equation (63),and the integration with an arbitrary adjoined weighting function ω^(†)yields the following expression:

$\begin{matrix}{{{\delta\mu}\lbrack{\delta\vartheta}\rbrack} \equiv {\frac{\langle\left. \omega^{\dagger} \middle| {c_{1}\varphi} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\vartheta}} \right.\rangle}{\delta\vartheta}}} & (65)\end{matrix}$

As the variation of the control parameter Θ is such that the secondneutron eigenvalue μ is driven towards 1, i.e. that μ+δμ is driventowards 1, the approximation δμ≡1−μ is imposed. With this approximationand the Equation (65), the following equation is obtained:

$\begin{matrix}{{\delta\vartheta} \cong {\frac{\langle\left. \omega^{\dagger} \middle| {c_{1}\varphi} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\vartheta}} \right.\rangle}\left( {1 - \mu} \right)}} & (66)\end{matrix}$

Using the equivalence μc₁φ=j_(out)−{circumflex over (Ω)}Ŷj_(out), δΘ isgiven by:

$\begin{matrix}{{{\delta\vartheta} \cong \frac{\langle\left. \omega^{\dagger} \middle| {j_{out} - {c_{1}\varphi} + {\hat{\Omega}\; \hat{Y}j_{out}}} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\vartheta}} \right.\rangle}} = \frac{\langle\left. \omega^{\dagger} \middle| r_{out} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\vartheta}} \right.\rangle}} & (67)\end{matrix}$

where r_(out) represents the residual of the interface outcurrentbalance equation, i.e. Equation (39), with the imposed target value forthe second neutron eigenvalue μ equal to 1.

The control parameter Θ is then updated such that:

$\begin{matrix}{\vartheta^{({new})} = {\vartheta^{({prev})} + \frac{\langle\left. \omega^{\dagger} \middle| r_{out}^{({prev})} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\vartheta}} \right.\rangle}}} & (68)\end{matrix}$

It is assumed that

$\frac{\partial\varphi}{\partial\vartheta} \cong \frac{\partial\overset{\sim}{\varphi}}{\partial\vartheta}$

which will be sufficient for a successful application of the variationalprinciple to the control parameter Θ.

For computing

$\frac{\partial\overset{\sim}{\varphi}}{\partial\vartheta},$

the above described conditioning with application of the Taylor formulais also used in the same manner.

The eighth aspect of the invention offers the same advantages as theabove described seventh aspect of the invention, and may be used incomplement of the first, second or third aspects of the invention in amanner similar to that for the above described seventh aspect of theinvention, in order to improve the computational efficiency.

FIG. 7 represents a set of convergence curves for different eigensystemiterative solving procedures with the number of pursued NEM iterativesteps in abscissa and the fast flux error in ordinate, which is themaximum solution error of the neutron flux for the fast neutrons energygroup. Since all energy groups are spectrally coupled, the maximumsolution errors of the neutron flux for the other energy groups arequite similar to the fast flux error. Thus, the comparison betweendifferent eigensystem iterative solving procedures is possible throughthese convergence curves.

The curve 50 is the convergence curve for the classic Coarse MeshRebalancing (CMR) procedure with use of single precision, and the curve51 is the convergence curve for the classic CMR procedure with use ofdouble precision. The single and double precisions are the classicprecisions used for representing the floating point numbers in acomputer program product.

The curve 52 is the convergence curve for the iterative procedureaccording to the first and the seventh aspects of the present inventionfor a single energy group with use of single precision, and the curve 53is the convergence curve for the iterative procedure according to thesame aspects of the present invention with use of double precision.

The curve 54 is the convergence curve for the iterative procedureaccording to the first and the second aspects of the present inventionfor four coarse energy groups with use of single precision, and thecurve 55 is the convergence curve for the iterative procedure accordingto the same aspects of the present invention with use of doubleprecision.

The curve 56 is the convergence curve for the iterative procedureaccording to the first, the second and seventh aspects of the presentinvention for four coarse energy groups with use of single precision,and the curve 57 is the convergence curve for the iterative procedureaccording to the same aspects of the present invention with use ofdouble precision.

The procedure according to the first, the second and seventh aspects forfour coarse energy groups with use of double precision (curve 57) is theprocedure with both the smallest error with a value substantially equalto 10⁻¹³ for 50 pursued NEM iteration steps and the best computationalefficiency, as illustrated by the gradient of the convergence curveswhich is maximum for the curve 57.

The procedure according to the first and the seventh aspects for asingle energy group with use of double precision (curve 53) alsoprovides a small error with a value substantially equal to 10⁻⁹ for 130pursued NEM iteration steps. Moreover, the computational efficiency ofthis procedure (curve 53) is worth than the computational efficiency ofthe procedure according to the first and the second aspects for fourcoarse energy groups with use of double precision (curve 55), becausethe gradient of the curve 55 is greater than the gradient of curve 53.

The procedure according to the first and the second aspects for fourcoarse energy groups with use of double precision provides an errorsubstantially equal to 10⁻⁷ for 40 pursued NEM iteration steps.

Thus, the double precision has a great effect on the value of the fastflux error, but has no impact on the computational efficiency of themodelizing method.

The seventh aspect of the present invention has also a greater effect onthe value of the fast flux error, than on the computational efficiency.

The first and the second aspects for four coarse energy groups providesthe best computational efficiency, and FIG. 7 illustrates that thegreatest gradient for an error value comprised between 10⁰ and 10⁻⁵ isobtained for the curves 54 to 57, which all correspond to the first andthe second aspects.

As set forth in the previous description, the first to eighth aspects ofthe invention help to achieve robust modelizing methods providing abetter computational efficiency so that relevant core simulations can beobtained within short computational time period.

It should be kept in mind that the first aspect in itself helps inachieving this goal and can thus be used separately from the second toeighth aspects. Also, the second, third, seventh and eighth aspects arenot necessarily implemented with the first aspect or with any one of thefourth to sixth aspects.

1-10. (canceled)
 11. A computer implemented method for modelizing anuclear reactor core, comprising the steps of: partitioning the core incubes to constitute nodes of a grid for computer implementedcalculation, calculating neutron flux by using an iterative solvingprocedure of at least one eigensystem, components of an iterant of theeigensystem corresponding either to a neutron flux, to a neutronoutcurrent or to a neutron incurrent, for a respective cube to becalculated, wherein a control parameter is varied to impact a neutroneigenvalue μ through a perturbed interface current equation and drivethe neutron eigenvalue μ towards
 1. 12. The method of claim 11 whereinthe control parameter Θ is given by:$\vartheta^{({new})} = {\vartheta^{({prev})} + \frac{\langle\left. \omega^{\dagger} \middle| r_{out}^{({prev})} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\vartheta}} \right.\rangle}}$wherein ω^(†) is a weighting function, Φ designates the neutron flux,r_(out) is given by: j_(out)−c₁φ+{circumflex over (Ω)}Ŷj_(out), Ŷ is acoupling operator, j_(out) designates the neutron outcurrent,{circumflex over (Ω)} is a mono-directional current throughflowoperator, and c₁ is an operator depending on the cube properties. 13.The method of claim 12 wherein the control parameter is a boronconcentration.
 14. The method of claim 12 wherein the control parameteris a parameter for controlling a rod insertion depth for a group ofselected control rods.
 15. The method of claim 11 wherein the controlparameter is a neutron eigenvalue λ, and wherein, during an iteration ofthe iterative solving procedure, the neutron eigenvalue λ is varied todrive the neutron eigenvalue μ towards 1, and wherein the variation ofthe neutron eigenvalue λ is an application of a variation δλ to λ suchthat λ is driven to λ+δλ, which impacts the neutron eigenvalue μ throughthe perturbed interface current equation.
 16. The method of claim 15,wherein the neutron eigenvalue λ is preconditioned prior to theiteration of the eigensystem solving procedure according to:$\lambda^{({new})} = {\lambda^{({prev})} + \frac{\langle\left. \omega^{\dagger} \middle| r_{out}^{({prev})} \right.\rangle}{\langle\left. \omega^{\dagger} \middle| {c_{1}\frac{\partial\varphi}{\partial\lambda}} \right.\rangle}}$where ω^(†) is a weighting function, Φ designates the neutron flux,r_(out) is given by: j_(out)−c₁φ+{circumflex over (Ω)}Ŷj_(out), Ŷ is acoupling operator, j_(out) designates the neutron outcurrent,{circumflex over (Ω)} is a mono-directional current throughflowoperator, and c₁ is an operator depending on the cube properties. 17.The method of claim 16 wherein the partial operator derivative$\frac{\partial\varphi}{\partial\lambda}$ is iteratively solved throughapplication of:$\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda} \right\rbrack^{({new})} = {s + {2\lambda \; {\hat{B}\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\vartheta} \right\rbrack}^{({old})}} - {\lambda^{2}{{\hat{B}}^{2}\left\lbrack \frac{\partial\overset{\sim}{\varphi}}{\partial\vartheta} \right\rbrack}^{({old})}}}$${{where}\mspace{14mu} \frac{\partial\varphi}{\partial\lambda}} \cong \frac{\partial\overset{\sim}{\varphi}}{\partial\lambda}${circumflex over (B)}=({circumflex over (R)}−Ŝ _(LD))⁻¹ {circumflex over(F)}s={circumflex over (B)}({circumflex over (R)}−Ŝ _(LD))⁻¹ ĜŶj_(out)=({circumflex over (R)}−Ŝ _(LD))⁻¹ {circumflex over(F)}({circumflex over (R)}−Ŝ _(LD))⁻¹ ĜŶj _(out) {circumflex over (R)}designates a removal operator, Ŝ_(LD) is a lower-diagonal part of aninscatter operator Ŝ, {circumflex over (F)} is a production operator,and Ĝ is an incurrent operator.
 18. The method of any one of claims 11further comprising a step of building a nuclear reactor core on thebasis of the calculated neutron flux.
 19. The method of claim 11 furthercomprising a step of operating a nuclear reactor core on the basis ofthe calculated neutron flux.
 20. A computer program product residing ona computer readable medium and comprising a computer program for runningon a computer the method according to claim 11.